library(Seurat)
library(dplyr)
library(patchwork)
library(tidyverse)
library(data.table)
library(readr)
library(harmony)
library(celldex)
library(SingleR)
library(BiocManager)
library(Matrix)
pathHC1<-"D:/R语言/脓毒数据HC1"
sepsisHC1<-Read10X(pathHC1)
sepsisHC1
pathHC2<-"D:/R语言/脓毒数据HC2"
sepsisHC2<-Read10X(pathHC2)
pathSep1<-"D:/R语言/脓毒数据Sep1"
sepsisSep1<-Read10X(pathSep1)
pathSep2<-"D:/R语言/脓毒数据Sep2"
sepsisSep2<-Read10X(pathSep2)
pathSep20<-"D:/R语言/脓毒数据Sep20"
sepsisSep20<-Read10X(pathSep20)
pathSep30<-"D:/R语言/脓毒数据Sep30"
sepsisSep30<-Read10X(pathSep30)
pathHC20<-"D:/R语言/脓毒数据HC20"
sepsisHC20<-Read10X(pathHC20)
pathHC30<-"D:/R语言/脓毒数据HC30"
sepsisHC30<-Read10X(pathHC30)
pathHC27<-"D:/R语言/脓毒数据HC27"
sepsisHC27<-Read10X(pathHC27)
gene_expressionHC1 <- sepsisHC1[["Gene Expression"]]
antibody_capture <- sepsisHC1[["Antibody Capture"]]
sepsisHC1<- CreateSeuratObject(counts = gene_expressionHC1, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsis[["ADT"]] <- CreateAssayObject(counts = antibody_capture)
sepsisHC1[["percent.mt"]] <- PercentageFeatureSet(sepsisHC1, pattern = "^MT-")
VlnPlot(sepsisHC1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
sepsisHC1 <- subset(sepsisHC1, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionHC2<-sepsisHC2[["Gene Expression"]]
sepsisHC2<-CreateSeuratObject(counts = gene_expressionHC2, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisHC2[["percent.mt"]] <- PercentageFeatureSet(sepsisHC2, pattern = "^MT-")
sepsisHC2 <- subset(sepsisHC2, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionSep1<-sepsisSep1[["Gene Expression"]]
sepsisSep1<-CreateSeuratObject(counts = gene_expressionSep1, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisSep1[["percent.mt"]] <- PercentageFeatureSet(sepsisSep1, pattern = "^MT-")
sepsisSep1 <- subset(sepsisSep1, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionSep2<-sepsisSep2[["Gene Expression"]]
sepsisSep2<-CreateSeuratObject(counts = gene_expressionSep2, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisSep2[["percent.mt"]] <- PercentageFeatureSet(sepsisSep2, pattern = "^MT-")
sepsisSep2 <- subset(sepsisSep2, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionSep20<-sepsisSep20[["Gene Expression"]]
sepsisSep20<-CreateSeuratObject(counts = gene_expressionSep20, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisSep20[["percent.mt"]] <- PercentageFeatureSet(sepsisSep20, pattern = "^MT-")
sepsisSep20 <- subset(sepsisSep20, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionSep30<-sepsisSep30[["Gene Expression"]]
sepsisSep30<-CreateSeuratObject(counts = gene_expressionSep30, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisSep30[["percent.mt"]] <- PercentageFeatureSet(sepsisSep30, pattern = "^MT-")
sepsisSep30 <- subset(sepsisSep30, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionHC20<-sepsisHC20[["Gene Expression"]]
sepsisHC20<-CreateSeuratObject(counts = gene_expressionHC20, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisHC20[["percent.mt"]] <- PercentageFeatureSet(sepsisHC20, pattern = "^MT-")
sepsisHC20 <- subset(sepsisHC20, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionHC27<-sepsisHC27[["Gene Expression"]]
sepsisHC27<-CreateSeuratObject(counts = gene_expressionHC27, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisHC27[["percent.mt"]] <- PercentageFeatureSet(sepsisHC27, pattern = "^MT-")
sepsisHC27 <- subset(sepsisHC27, subset = nFeature_RNA > 200 & percent.mt < 10)
gene_expressionHC30<-sepsisHC30[["Gene Expression"]]
sepsisHC30<-CreateSeuratObject(counts = gene_expressionHC30, project = "pbmc3k", min.cells = 3, min.features = 200)
sepsisHC30[["percent.mt"]] <- PercentageFeatureSet(sepsisHC30, pattern = "^MT-")
sepsisHC30 <- subset(sepsisHC30, subset = nFeature_RNA > 200 & percent.mt < 10)
samlist<-list(sepsisHC1, sepsisHC2, sepsisHC20, sepsisHC27, sepsisHC30, sepsisSep1, sepsisSep2, sepsisSep20, sepsisSep30)
sepsisHC1$sample <- "HC1"
sepsisHC2$sample <- "HC2"
sepsisHC20$sample <- "HC20"
sepsisHC27$sample <- "HC27"
sepsisHC30$sample <- "HC30"
sepsisSep1$sample <- "Sep1"
sepsisSep2$sample <- "Sep2"
sepsisSep20$sample <- "Sep20"
sepsisSep30$sample <- "Sep30"


sepsis <- merge(samlist[[1]], y = samlist[-1], add.cell.ids = c("HC1","HC2","HC20","HC27","HC30","Sep1","Sep2","Sep20","Sep30"))
all_cells <- rownames(sepsis@meta.data)
sample_labels <- sapply(all_cells, function(cell) {
  prefix <- gsub("_.*", "", cell)
  switch(prefix,
         "HC1" = "HC1",
         "HC2" = "HC2",
         "HC20" = "HC20",
         "HC27" = "HC27",
         "HC30" = "HC30",
         "Sep1" = "Sep1",
         "Sep2" = "Sep2",
         "Sep20" = "Sep20",
         "Sep30" = "Sep30",
         "Unknown")
})
sepsis$sample <- sample_labels
sepsis <- NormalizeData(sepsis, normalization.method = "LogNormalize", scale.factor = 10000)
sepsis <- FindVariableFeatures(sepsis, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(sepsis), 10)
top10
plot1 <- VariableFeaturePlot(sepsis)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
all.genes <- rownames(sepsis)
sepsis <- ScaleData(sepsis) 
head(sepsis@meta.data) 
unique(sepsis@meta.data$orig.ident)
sepsis <- RunPCA(sepsis, features = VariableFeatures(object = sepsis))
DimPlot(sepsis, reduction = "pca") 
ElbowPlot(sepsis, ndims = 50)
sepsis <- FindNeighbors(sepsis, dims = 1:50)
sepsis <- FindClusters(sepsis, resolution = 0.5)
head(Idents(sepsis), 5)
sepsis <- RunUMAP(sepsis, group.by = "sample", dims = 1:50)
DimPlot(sepsis, reduction = "umap", group.by = "sample")

#手动标记
sepsis <- JoinLayers(sepsis)
sepsis.markers <- FindAllMarkers(sepsis, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
head(sepsis.markers, n=5)
top_markers <- sepsis.markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC)
print(top_markers, n=50)
mean(GetAssayData(sepsis, assay = "RNA", layer = "counts")["MS4A6A", sepsis$seurat_clusters == "0"] > 0)
VlnPlot(sepsis, features = c("MS4A6A", "LGALS2", "TMEM176B", "CD93", "LYZ"))
FeaturePlot(sepsis, features = c("MS4A6A", "LGALS2", "TMEM176B", "CD93", "LYZ"))
for (cl in unique(sepsis$seurat_clusters)) {
  cat("\n=== Cluster", cl, "Top 5 markers ===\n")
  print(head(subset(sepsis.markers, cluster == cl), 5))
}
#显示每个cluster的高表达基因
cluster_annot <- c(
  "0"  = "Inflammatory Monocyte",
  "1"  = "Neutrophil",
  "2"  = "Naive CD4+ T cell",
  "3"  = "Activated Monocyte/DC",
  "4"  = "Non-classical Monocyte",
  "5"  = "Megakaryocyte/Platelet",
  "6"  = "NK cell",
  "7"  = "Cytotoxic NK cell",
  "8"  = "CD8+ Effector T cell",
  "9"  = "Neutrophil (Granule-rich)",
  "10" = "Inflammatory Monocyte",
  "11" = "NK cell",
  "12" = "cDC2",
  "13" = "Naive CD8+ T cell",
  "14" = "B cell",
  "15" = "Erythroid Progenitor",
  "16" = "Erythrocyte",
  "17" = "Plasma cell",
  "18" = "pDC",
  "19" = "Low Quality / Stress",
  "20" = "Unknown / Non-immune",
  "21" = "Proliferating cell",
  "22" = "cDC1"
)
sepsis <- RenameIdents(sepsis, cluster_annot)
DimPlot(sepsis, label = TRUE, repel = TRUE) + NoLegend()
VlnPlot(sepsis, features = c("CD3D", "CD4", "CD8A", "CD14", "CD19", "CD20"), ncol = 4)
top10 <- sepsis.markers %>%
  group_by(cluster) %>%
  filter(avg_log2FC > 1) %>%
  slice_head(n = 10) %>%
  ungroup()
DoHeatmap(sepsis, features = top10$gene) + 
  NoLegend() +
  theme(
    axis.text.x = element_text(size = 6, angle = 60, hjust = 1),
    axis.text.y = element_text(size = 7),
    plot.margin = unit(c(0.5, 0.5, 1.5, 0.5), "cm")
  )

#自动标记
hpca.se <- HumanPrimaryCellAtlasData()
sce <- as.SingleCellExperiment(sepsis)
pred.hpca <- SingleR(
  test = sce,
  ref = hpca.se,
  labels = hpca.se$label.main,
  assay.type.test = "logcounts"  # 注意：Seurat 默认存储在 "data" slot，对应 log-normalized
)
table(pred.hpca$labels)
plotScoreHeatmap(pred.hpca)
